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Abstract 

This  paper  addresses  the  development  of  parameter  estimation  techniques  for  a  class  of  models 
used  to  characterize  hysteresis  and  constitutive  nonlinearities  inherent  to  ferroelectric,  ferromagnetic 
and  ferroelastic  materials  employed  in  a  wide  range  of  actuators  and  sensors.  These  models  are 
formulated  as  integral  equations  with  known  kernels  and  unknown  densities  to  be  identified  through 
least  squares  techniques.  Due  to  the  compactness  of  the  integral  operators,  the  resulting  discretized 
models  inherit  ill-posedness  which  must  be  accommodated  through  regularization.  The  accuracy  of 
regularized  finite-dimensional  models  is  illustrated  through  comparison  with  experimental  data. 
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1  Introduction 


Hysteresis  and  constitutive  nonlinearities  are  inherent  properties  of  a  wide  range  of  piezoceranric 
(PZT),  magnetic,  and  shape  memory  alloy  (SMA)  compounds  being  considered  as  actuators  and  sen¬ 
sors  in  aeronautic,  aerospace,  automotive,  industrial  and  biomedical  applications.  In  some  regimes, 
hysteresis  can  be  mitigated  through  restricted  input  levels,  amplifier  designs,  or  feedback  mech¬ 
anisms.  In  general,  however,  these  properties  are  ubiquitous  and  are  often  inexorably  related  to 
material  attributes  which  provide  the  materials  with  unique  transducer  capabilities.  For  example, 
damping  provided  by  a  shape  memory  tendon  is  proportional  to  the  area  of  the  hysteresis  loop. 
Hence  optimal  tendon  design  to  attenuate  wind  or  earthquake-induced  vibrations  in  civil  structures 
requires  maximal  hysteresis  rather  than  operation  in  approximately  linear  regimes.  This  necessi¬ 
tates  the  development  of  models,  parameter  estimation  techniques,  and  model-based  control  designs 
which  incorporate  the  nonlinear  hysteresis  mechanisms  in  a  manner  which  facilitates  subsequent 
implement  ation . 

A  number  of  modeling  strategies  have  been  proposed  but  three  stand  out  in  the  sense  that  they 
provide  unified  frameworks  for  characterizing  hysteresis  in  ferroelectric,  ferromagnetic  and  ferroelas- 
tic  materials,  which  are  collectively  referred  to  as  ferroic  compounds.  These  three  approaches  are 
homogenized  energy  models  [8,10,18-20,24,25],  Preisach  formulations  [2,3,5,7,16,26],  and  domain 
wall  models  [6,8,14,21,22], 

As  detailed  in  [17],  the  domain  wall  framework  is  efficient  to  implement  but  requires  a  priori 
knowledge  of  turning  points  to  guarantee  closure  of  biased  minor  loops.  This  precludes  its  use  in 
certain  feedback  designs  where  turning  points  are  determined  by  measured  or  estimated  states  of  the 
system.  The  classical  Preisach  framework  requires  the  properties  of  congruency  and  deletion  [11, 12] 
which  are  overly  restrictive  for  a  number  of  materials  and  applications.  As  detailed  in  [1],  these 
restrictions  have  been  eliminated  for  magnetic  material  characterization  through  the  development  of 
extended  Preisach  models;  however,  these  extensions  come  at  the  price  of  increased  complexity,  and 
the  maturity  of  extended  Preisach  models  for  PZT  and  SMA  lags  far  behind  the  magnetic  theory. 
The  homogenized  energy  framework  is  the  most  recent  of  the  three  and  is  based  on  a  combina¬ 
tion  of  energy  analysis  at  the  lattice  level  and  stochastic  homogenization  techniques  to  construct 
macroscopic  models.  Due  to  its  energy  basis,  it  incorporates  certain  frequency,  temperature,  and 
stress  dependencies  which  makes  it  applicable  for  a  wide  range  of  transducer  designs  and  operating 
regimes.  Furthermore,  it  is  illustrated  in  [23]  that  this  framework  provides  an  energy  basis  for  certain 
extended  Preisach  formulations.  Due  to  its  flexibility  and  generality,  we  employ  the  homogenized 
energy  framework  in  this  paper  and  construct  parameter  estimation  techniques  in  this  context. 

In  Section  2,  we  summarize  the  homogenized  energy  model  for  ferroelectric  (e.g.,  PZT),  ferro¬ 
magnetic  (e.g.,  iron,  Terfenol-D),  and  ferroelastic  (e.g.,  SMA)  compounds  to  illustrate  its  structure 
and  generality.  Specifically,  it  will  be  noted  that  the  models  are  formulated  as  integral  equations 
with  known  kernels  and  unknown  densities.  The  compactness  of  the  integral  operator  is  established 
in  Section  3  and,  in  Section  4,  it  is  shown  that  this  leads  to  ill-posedness  in  the  inverse  problem 
associated  with  estimating  the  densities  given  a  set  of  measured  input-output  data.  Experimental 
validation  results  are  summarized  in  Section  5.  To  simplify  the  discussion,  Sections  3  through  5  will 
focus  on  the  polarization  model.  Due  to  the  general  nature  of  the  framework,  however,  the  results 
are  equally  applicable  to  the  magnetic  and  shape  memory  alloy  models. 


1 


2  Free  Energy  Framework 

We  summarize  in  this  section  the  free  energy  framework  for  characterizing  hysteresis  in  ferroelectric, 
ferromagnetic  and  ferroelastic  materials.  Details  regarding  the  development  and  validation  of  these 
models  can  be  found  in  [10, 19,20,25]. 

2.1  Polarization  Model 

2.1.1  Local  Polarization  Model 

To  quantify  the  internal  and  electrostatic  energy  at  the  lattice  level,  we  employ  the  Helmholtz  energy 

'  ^(P  +  Pr)2  ,  P  <  —Pi 

if>(P)  =  <  h(p  -  pR )2  ,  P>  Pi  (1) 

k  Hpi  ~  Pn)  (f  -  PR )  t  1^1  <  pi- 

and  Gibbs  energy 

G(E,P)=^(P)-EP  (2) 

where  E  and  P  are  the  electric  field  and  polarization  and  Pj  and  Pr  denote  the  inflection  point  and 
polarization  at  which  the  minimum  of  if)  occurs. 

In  the  absence  of  thermal  activation,  minimization  of  G  for  fixed  field  inputs  yields  the  local 
polarization  relation 

P(E)  =  E+  SPR  (3) 

V 

where  the  parameter  5  has  a  value  of  1  for  positively  oriented  dipoles  and  —1  for  negative  orientations. 
To  specify  5  and  hence  P  more  specifically  in  terms  of  the  initial  dipole  orientations  and  previous 
switches,  we  employ  Preisach  notation  and  take 

'  [P(£;Pc,<£)](0)  ,  r(t)  =  0 

[P(P;PC,  £)](£)  =  <  p  -  Pr  ,  r(t)  /  0  and  P(maxr(f))  =  -Ec  (4) 

,  f  +Pr  ,  r{t)  /  0  and  £(maxr(t))  =  Ec. 

Here 

'  f  -pR  ,  E( 0)  <  -Ec 

[P(E;Ec,0m  =  C  ,  -Ec  <  E{ 0)  <  Ec 

.  f  +  pR  ,  E( 0)  >  Ec 

denotes  the  initial  dipole  distribution  and  transition  times  are  designated  by 

r(t)  =  {t  £  (0,  Tf]  |  E(t)  =  -Ec  or  E(t)  =  Ec} 
where  Tf  denotes  the  final  time  under  consideration.  The  local  coercive  field 

Ec  =  v(pR  ~  Pi)  (5) 

quantifies  the  field  at  which  the  negative  well  ceases  to  exist  and  hence  a  dipole  switch  occurs.  As 
shown  in  Figure  1,  the  kernel  (3)  or  (4)  quantifies  the  abrupt  hysteresis  transitions  associated  with 
homogeneous,  single  crystal  compounds  in  the  absence  of  thermal  activation. 
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Figure  1:  Local  polarization  P  given  by  equation  (6)  with  high  thermal  activation  ( - )  and  limiting 

polarization  P  specified  by  (4)  in  the  absence  of  thermal  activation  ( - ). 


To  incorporate  thermal  activation,  the  Gibbs  energy  G  and  relative  thermal  energy  kT/V  are 
balanced  through  the  Boltzmann  relation 

H(G)  =  Ce-GV'kT. 


Here  k  is  Boltzmann’s  constant,  V  is  a  reference  volume,  T  is  the  temperature  in  degrees  Kelvin, 
and  C  is  a  constant  chosen  to  ensure  that  integration  over  all  admissible  inputs  yields  a  probability 
fi  of  unity.  As  detailed  in  [25],  the  kernel  in  this  case  is 

P  =  x+  (P+)  +  x-  ( P- )  (6) 


where  x+  and  x_  denote  the  fractions  of  dipoles  having  positive  and  negative  orientations  and  (P+) 
and  (P~)  denote  the  average  expected  polarizations  associated  with  the  two  orientations.  Since  the 
expected  polarization  values  are  obtained  by  integrating  the  product  P/j,(G(P))  over  all  admissible 
configurations,  it  follows  that 


(P+) 


Pe-G(E,P,T)V/kTdp 


/  e-G(E,P,T)V/kTdp 

I  Pi 


r-Pi 


Pe-G(E,P)V/kTdp 


(P-)  = 


r-Pi 


-G(E,P)V/kTdp 


The  denominator  results  from  the  evaluation  of  the  integration  constant  C  whereas  it  is  illustrated 
in  [25]  that  the  use  of  the  inflection  points  =t Pj,  to  simplify  evaluation  of  the  integrals,  can  be  justified 
though  either  asymptotic  analysis  or  energy  arguments. 

Debye  arguments  yield  the  differential  equations 


x+  =  -p+-X+  +p-+x _ 

X-  =  —p _ |_a:_  +  p. | _ x+ 

quantifying  the  evolution  of  the  respective  dipole  fractions.  For  implementation,  these  relations  can 
be  simplified  to  the  single  differential  equation 


X+  =  -P+-X+  +  p-  +  (l  -  x+) 


through  the  identity  x+  +  X-  =  1. 
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The  likelihoods  of  switching  from  positive  to  negative  orientations,  or  vice  versa,  are  respectively 
quantified  by 


rPi+e 


e-G(E,PI,T)V/kTdp 


P+-  = 


I  Pi 


T(T) 


-G(E,P,T)V/kTdp 


I  Pi 


P-+ 


1 

nf) 


G(E,-PI,T)V/kTdp 


G(E,P,T)V/kT  dp 


where  e  is  taken  to  be  a  small  positive  constant.  The  relaxation  term  T  quantifies  the  frequency 
at  which  jumps  are  attempted  whereas  the  remainder  of  the  definition  characterizes  the  probability 
of  achieving  the  energy  required  to  exit  respective  potential  wells.  It  is  detailed  in  [25]  that  this 
probability  increases  when  the  relative  thermal  energy  kT  /V  approaches  the  Gibbs  energy  G. 


Remark  1.  It  is  proven  in  [17,  25]  that  the  kernel  P  given  by  (6)  converges  to  the  piecewise  linear 
kernel  (3)  or  (4)  in  the  limit  kT /V  — >  0  of  negligible  thermal  activation.  It  also  follows  immediately 
that  P  given  by  (6)  satisfies 


P  <  P T 

—  mm 


p  >  p~ . 

—  mm 


- b  Pr  ,  positively  oriented  dipoles 

V 


- Pr  ,  negatively  oriented  dipoles 

V 


as  depicted  in  Figure  1.  For  all  fields  E  £  C[a,b]  with  a,b  finite,  it  thus  follows  that  P  £  L1(a,  b) 
and  P  £  L2(a,b)  for  P  given  by  (3),  (4)  or  (6). 


2.1.2  Macroscopic  Polarization  Model 

To  incorporate  the  effects  of  material  nonhomogeneities,  poly  crystallinity,  and  variable  effective  fields 
Ee  =  E  +  Ej ,  where  Ej  denotes  interaction  fields,  we  assume  that  Ej  and  the  local  coercive  field  Ec, 
defined  in  (5),  are  manifestations  of  underlying  distributions  rather  than  constants.  This  yields  the 
macroscopic  model 


ip(mt) 


r*oo  /*oo 


'  0  J  —  oo 


MEMErfiP^E  +  En  Ec,0}(t)  dEj  dEc 


r*oo  /*oo 


'  0  J  —  oo 


v[Ec,  Ej)[P{E  +  Ei;  Ec,m)  dE 7  dEc 


(7) 


where  v\  and  are  densities  associated  with  Ec  and  Ej. 

To  accommodate  physical  criteria,  we  assume  that  >  0  and  z^2  >  0  satisfy  the  conditions 


(i)  u\  (x)  defined  for  x  >  0, 

(ii)  P2(-x)  =  iy2(x), 

(iii)  Mx)|  <  Cle~aix, 
h(x)j  <  c2e~a2^ 


for  positive  ci,ai,C2,a2.  The  restricted  domain  in  (i)  reflects  the  fact  that  the  coercive  field  Ec  is 
positive  whereas  the  symmetry  enforced  in  the  interaction  field  through  (ii)  yields  the  symmetry 
observed  in  low- field  Rayleigh  loops.  Hypothesis  (iii)  incorporates  the  physical  observation  that  the 
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coercive  and  interaction  fields  decay  as  a  function  of  distance  and  guarantees  that  integration  against 
the  piecewise  linear  kernel  yields  finite  polarization  values. 

Remark  2.  It  is  illustrated  in  [17,20]  that  for  certain  operating  regimes,  reasonable  accuracy  can  be 
obtained  using  a  priori  choices  for  and  1*2-  Motivated  by  densities  employed  in  Preisach  models 
for  magnetic  compounds,  one  such  choice  is 

ME,)  =  Cie-MEc/Ec)/2c]* 

MEi )  =  c2e_E?/262. 

For  general  operating  regimes,  however,  identification  of  general  densities  u\  and  V2  is  necessary  to 
achieve  accurate  material  characterization  throughout  both  biased  and  unbiased  operating  regimes. 

Remark  3.  Formulation  of  the  model  in  terms  of  the  joint  density  EC,E j)  =  ui(Ec)i>2(Ei) 
increases  the  dimensionality  of  the  parameter  estimation  problem  but  yields  models  that  exhibit  a 
linear  dependence  on  parameters.  This  proves  advantageous  when  constructing  linear  adaptive  control 
techniques  and  for  this  reason,  we  focus  heavily  on  this  case  in  subsequent  analysis. 


(9) 


2.1.3  Discretized  Polarization  Model 

For  implementation  purposes,  Gaussian  or  Newton-Cotes  quadrature  routines  can  be  employed  to 
approximate  the  integrals,  thus  yielding  the  system 

Ni  Nj 

[p(E)m  =  +E-,Eci,omviwj  ao) 

i= 1  3= 1 

where  Vj  and  Wj  are  the  weights  associated  with  the  quadrature  rules  and  ECi,  Ej.  are  the  abscissas. 
Highly  efficient  algorithms  for  implementing  the  approximate  model  (10)  with  general  densities  v\ 
and  V2  are  provided  in  [17,25]. 

Formulation  of  the  model  in  terms  of  the  joint  density  v  =  v\  ■  zz2  increases  the  dimensionality  of 
the  parameter  estimation  problem  from  Nt  +  Nj  to  Nt  ■  Nj  but  yields  a  system  which  depends  linearly 
on  the  parameters  {n(ECi,  Ejj)}.  This  permits  implementation  of  linear  least  squares  algorithms  and 
linear  adaptive  identification  and  control  techniques. 

To  formulate  the  discretized  model  (10)  as  a  linear  system  in  terms  of  v  =  u\  ■  1/2,  we  define  the 
Ni  x  Nj  matrices  A(E)  and  $  to  have  components 

[A(E)]ij=  - -Jl  +  PR5(E;ECi,EM  viWj 

[Q\ij  =  v(ECi ,  Ejj ) . 

For  N  =  Ni  ■  Nj,  we  define  the  JV  x  1  vector  q  and  1  x  N  vector  a(E)  by 

q  =  vec  (Q)  ,  a(E)  =  [vec  (A(E))]T 

where  ‘vec’  denotes  the  vector  concatenation  of  the  respective  matrices.  The  discretized  polarization 
model  (10)  can  then  be  formulated  as  the  linear  system 

P(E)  =  a(E)q.  (11) 

We  note  that  r/  is  considered  known  and  fixed  in  this  formulation  and  is  incorporated  in  a(E). 
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2.2  Magnetization  Model 

The  magnetization  model  is  analogous  so  we  simply  summarize  it  here.  Details  regarding  its  devel¬ 
opment  can  be  found  in  [18, 19]. 

The  hysteretic  relation  between  the  magnetic  field  H  and  magnetization  M  in  ferromagnetic 
materials  is  characterized  by  the  model 

roc  roc 

[M(H)}(t)=  /  /  v\(H(?)v2{Hi)\M(H  +  Hj;  Hc,  £)](£)  dHj  dHc.  (12) 

J  0  J —oo 


In  the  absence  of  relaxation  processes,  the  kernel  M  has  the  form 


M  = 


- b  5  Mr 

V 


whereas  it  is  given  by 


M  =  x+  (M+)  +  x _  (M_) 


when  the  Gibbs  and  relative  thermal  energies  are  balanced  to  incorporate  thermal  relaxation  pro¬ 
cesses.  The  dipole  fractions  x+,x_,  average  magnetizations  (M+) ,  (M_),  local  remanence  magne¬ 
tization  Mr,  and  switching  parameter  5  are  defined  in  a  manner  analogous  to  their  polarization 
counterparts  as  detailed  in  Section  2.1.  As  with  the  polarization  model,  the  goal  in  the  parameter 
identification  problem  is  to  estimate  v  given  data  measurements  {H^,  M *,},  k  =  1, . . . ,  Nj. 


2.3  Shape  Memory  Alloy  Model 

The  relation  between  input  stresses  a  and  strains  e  generated  in  shape  memory  alloy  compounds  as 
they  undergo  martensitic  phase  transformations  is  quantified  by  the  integral  relation 

roc  roc 

[s(a,T)}(t)  =  /  /  i^i{crR)i^2(crI)[£(a  +  07,  T;  ctr,  £)](t)  doj  daR  (13) 

J  0  J — 00 

which  has  the  same  general  form  as  the  polarization  model  (7)  and  magnetization  model  (12).  For 
regimes  in  which  thermal  activation  or  relaxation  mechanisms  are  significant,  the  kernel  e  is  given 
by 

e  =  x-  {£-)  +  x+  (e+)  +  xA  {sa) 

where  x_,x+  and  xa  respectively  denote  the  volume  fraction  of  martensite  minus,  martensite  plus, 
and  austenite  layers  in  1-D  SMA  compounds  and  (e_) ,  (e+)  and  (sa)  are  the  average  strains  asso¬ 
ciated  with  layers.  Details  regarding  the  construction  of  thermally  active  and  inactive  SMA  models 
can  be  found  in  [9, 10, 15, 17]. 


3  Compactness  of  the  Polarization  Operator 

The  behavior  of  the  inverse  problems  associated  with  estimating  the  densities  zq  and  v 2,  or  the  joint 
density  v  =  v\  ■  1^2,  is  intimately  related  to  properties  of  the  integral  operators  resulting  from  the 
models  (7),  (12)  and  (13).  In  this  section,  we  establish  that  the  operators  are  compact;  in  Section  4, 
we  use  this  fact  to  establish  the  ill-posedness  of  the  inverse  problem  associated  with  parameter 
estimation.  To  simplify  the  discussion,  we  focus  solely  on  the  polarization  model  in  subsequent 
sections.  Due  to  the  unified  nature  of  the  characterization  framework,  however,  the  results  also 
apply  to  the  magnetization  and  SMA  models. 
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We  focus  on  the  case  v  =  u\  ■  z^2,  which  yields  a  linearly  parameterized  problem  but  comes 
at  the  cost  of  increased  dimensionality.  The  kernel  P  in  (7)  can  be  taken  as  either  the  piecewise 
linear  relation  (3)  or  (4),  that  arises  when  thermal  activation  is  negligible,  or  the  kernel  (6)  which 
incorporates  thermal  processes.  The  important  observation  in  both  cases  is  the  property  that  for 
finite  input  fields  E  £  C[a,  b],  the  kernel  satisfies  P  £  L1(a,  b)  and  P  £  L2(a ,  b )  as  noted  in  Remark  1. 

It  was  noted  in  Section  1  that  the  homogenized  energy  framework  provides  an  energy  basis  for 
certain  extended  Preisach  models  and  hence  one  would  expect  similar  compactness  results  for  the 
two  theories.  The  reader  is  referred  to  Iyer  and  Shirley  [7]  for  theory  establishing  the  compactness 
of  the  classical  Preisach  operator. 

By  invoking  the  physical  decay  criteria  (8),  the  polarization  model  (7)  can  be  approximated  to 
arbitrary  accuracy  by  the  relation 

[P(E)](t)  =  ff  v(Ec,EI)\P(E  +  EI;Ec,S)](t)dEIdEc  (14) 


on  the  compact  domain 

=  {(Ec,  Ej )  £  M_|_  x  M  |  is(ECl  Ej)  P  e}. 

Furthermore,  we  let  the  minimum  and  maximum  admissible  input  fields  be  denoted  Emtn  and 
Emax  and  define 

-  -  [  —  [Emini  Emax\- 

We  consider  parameters  q  =  v  in  the  parameter  space 

Q  =  L2(n2)  (15) 

and  define  the  observation  operator  CP  =  P(E)  on  the  observation  space 

y  =  L2(Pmin,Pmax)-  (16) 

The  polarization  model  (7)  can  then  be  formulated  as 

y(E)  =  JCq(E) 


where 

E  £  Cpli]  C  L2(n i) 

and  the  parameter-to-observation  operator  K,  is  defined  by 

Kq=cJJk(.  +  EI,EMEc,E,)dEldEc. 
n2 


(17) 


From  the  property  that  P  given  by  (4)  or  (6)  satisfies  P  £  L1(Hi),  P  £  L2(Cl i),  as  noted  in 
Remark  1,  it  follows  that  k  £  L1(H)  and  k  £  L2(Q)  where 


0  —  X  ^2* 


The  property  that  k  £  T1(H)  is  typical  for  convolution  operators  whereas  k  £  L2(Q)  facilitates 
construction  of  a  generalized  Fourier  basis  for  the  operator.  We  employ  this  latter  property  to 
establish  that  K,  is  compact  for  the  given  choice  of  spaces.  As  a  prelude,  we  state  the  following 
theorem  which  is  Theorem  5.24.8  from  [13]. 


7 


Theorem  1.  Let  X  and  Y  be  Banach  spaces  and  let  /Cat  :  X  — >  Y,N  =  1,2, ,  be  a  sequence  of 
compact  linear  operators  converging  to  a  bounded  linear  operator  K.  :  X  —>  Y ;  that  is,  ||/C/v  —  /C||  — >  0 
as  N  — *  oo.  Then  1C  is  a  compact  linear  operator. 

Remark  4.  Consider  the  parameter  space  Q  and  observation  space  y  defined  in  (15)  and  (16).  The 
integral  operator  given  by  (17)  is  then  a  compact  operator.  We  establish  this  by  demonstrating  that 
K,  is  the  limit  of  a  sequence  of  finite  rank  operators  followed  by  the  use  of  Theorem  1. 


We  first  construct  an  orthonormal  basis  {fii}  for  L2(Ll).  It  is  illustrated  in  [13]  that 


Ms)  =  ,  „  =  exP 

V  f^max  -L^min 


o  •  n  S  Emin 

2m£- l-E 

u mm 


1  =  0,±1,±2,  ••• 


forms  an  orthonormal  basis  for  L2(fili).  With  an  analogous  basis  definition  for  L2(fl2)>  it  follows 
that  an  orthonormal  basis  for  L2(Q)  is 


fam(s,t,v)  =  Pi(s)ipm(t,v) 


which  we  re-index  as  {fii}. 

It  follows  that  every  /  G  L2(Q)  has  the  generalized  Fourier  series  representation 

/  =  E  fa)  fa 

i 

where  (■,  •}  denotes  the  usual  L2  inner  product.  The  norm  representation 


ii/ii2  =  Ei  </•*> 

i 


follows  from  Plancheral’s  theorem.  Moreover,  we  can  represent  K,  and  approximating  finite-rank 
operators  JCn  by 

£/  =  E  fa)  fa 

i 

N 

ICNf  =  X]  (/>  fa)  fa 

i=  1 


where  i pi  =  IC(j>i. 

To  establish  the  convergence  K,  — ►  /Cat,  we  note  that 


||/C/ -/Cat/H 


< 


< 


< 


E  (f>  fa)  & 

i>N+l 

E  \(f’fa)\\\fa 


i>N+ 1 


E 

^  i>iV+l 


(f,  fa) 


2  E 


L  i>7V+l 


E 


L  i>N+l 


where  the  third  inequality  follows  from  the  Schwartz  inequality.  Furthermore,  we  observe  that 


£ll*ll2 


,  L  Jn  i 


£ 

i 

y. 

i 

[ 

j  fji 


\Kct>i{E)\2  dE 


k(E  +  £7,  Ec)(pi(Ec,  E '/)  dEi  dE, 

'y  "  J  J  k(E  +  Ej ,  Ec)(f>i(Ec,  Ej)  dEj  dE, 
i  n2 

jf  |fc(£  +  £/,£c)|2d£/d£ 


dE 


dE 


dE  <  oo. 


^2 


The  last  step  follows  from  Plancheral’s  theorem.  Convergence  of  the  series  1 1 -0*  1 1 2  implies  that 

'y2i>N+ 1  llV'ill2  - ^  0  as  N  — ►  oo.  Thus  for  e  >  0,  there  exists  Ne  such  that  for  N  >  N£, 


||  /C  —  /Civil  =  sup 
/#  o 


\\JCf-JCNf\\ 

ll/ll 


<  e 


which  establishes  that 

lim  || K,  —  /Cjv ||  =  0. 

N^-oo 

Since  the  range  of  /Cat  is  finite,  it  follows  that  /Cat  is  a  compact  operator.  The  compactness  of  K. 
follows  from  Theorem  1  since  it  is  the  norm  limit  of  a  sequence  of  compact  operators. 

4  Parameter  Identification  Problem 

We  focus  on  the  problem  of  identifying  the  joint  density  v  in  the  polarization  model  (7)  or  density 
values  {u(ECi,  E /.)}  in  the  discretized  problem  (10)  due  to  the  fact  that  these  formulations  yield  the 
linear  parameterization  necessary  for  linear  adaptive  identification  or  control  techniques.  The  formu¬ 
lation  of  parameter  estimation  problems  associated  with  the  identification  of  independent  densities  vi 
and  1/2  or  parameters  arising  in  functional  representations  for  u\  and  v 2  —  e.g.,  C  =  c\  ■  C2,  E c,  c,  b  in 
(9)  —  is  analogous  but  requires  nonlinear  parameterizations  in  the  operator  and  matrix  formulations. 

For  the  operator  /C  defined  in  (17),  data  P  corresponding  to  input  fields  E  G  L2(Emin,  Emax), 
and  parameter  space  Q  =  L2(Q 2),  the  parameter  estimation  problem  can  be  formulated  as  follows: 
find  q  G  Q  so  that 

JCq  =  P.  (18) 

We  note  that  (18)  has  a  classical  solution  if  and  only  if  P  G  7 Z(JC),  where  7 Z(JC)  denotes  the  range 
of  1C,  which,  in  general,  will  not  be  true.  Instead  it  is  more  reasonable  to  consider  the  least  squares 
problem 

mmT(q),T(q)  =  h\ICq-Pfy.  (19) 

q&Q  l 

However,  because  /C  is  compact  with  infinite  dimensional  range,  the  Moore-Penrose  inverse  K)  is 
discontinuous  so  that  even  (19)  is  ill-posed  —  see  [4],  This  motivates  consideration  of  the  augmented 
functional 

Ta(q)=  l-\\lCq-P\\2y  +  aJ{q)  (20) 
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and  the  regularized  least  squares  minimization  problem 


min  Ta(q). 
q&Q 


(21) 


The  regularization  parameter  a  >  0  controls  the  tradeoff  between  goodness  of  fit  to  the  data  and 
stability  whereas  the  penalty  functional  J  provides  stability  and  allows  the  inclusion  of  a  priori 
information  regarding  the  parameter  q.  One  choice  for  J  is  the  Tikhonov  functional  which  we 
illustrate  in  the  context  of  the  discretized  problem. 

To  formulate  the  finite-dimensional  parameter  estimation  problem,  we  modify  the  linearly  pa¬ 
rameterized  system  (11)  to  reflect  measured  data.  We  define  the  TVj  x  Nj  matrices 


[Ak]ij 


Ek  +  Eij 
V 


+  PRS(Ek-  ECi,E j. 


ViWj 


[Q]ij  ~  v(ECi,  E^) 

and  vector  concatenations 

q  =  vec  (Q)  ,  ak  =  [vec(Afc)]T 

so  that  q  and  ak  are  respectively  1  x  IV  and  N  x  1  where  N  =  N 
vectors  V  and  V  are  defined  componentwise  by 

[P]k  =  P(Ek',  q)  ,  [P]k  =  Pk 

and  the  Nj  x  N  matrix  A  is  defined  row-wise  by 

[A]k  =  ak • 

The  discretized  polarization  model  (7)  can  then  be  formulated  as  the  linearly  parameterized  system 

V(Ek)  =  Aq 

—  see  [17,25]  for  highly  efficient  implementation  algorithms. 

The  unregularized  least  squares  problem  used  to  estimate  q  =  v  £  Q  =  WNi'Ni  given  measure¬ 
ments  {{Ek.  Pk)},  k  =  1, . . . ,  Nj  is  the  following: 

min T(q)  ,  T(q)  =  \  \\Aq-V\\2 

q&Q  2  (23) 

subject  to  qi  >  0,  j  =  1, . . . ,  N. 

Here  ||  •  ||  denotes  the  Euclidean  norm  in  To  incorporate  Tikhonov  regularization,  we  consider 
the  minimization  problem 

min Ta(q)  ,  T(q)  =  ]-\\Aq-Vk\\2  +  ^\\q\\2 

q&Q  2  2  (24) 

subject  to  qi  >  0,  j  =  1, . . . ,  N. 

Techniques  for  choosing  a  to  avoid  oversnroothing  solutions  as  well  as  a  solution  algorithm  for  (24) 
can  be  found  in  Vogel  [27] . 

We  do  not  consider  the  convergence  of  approximate  parameters  as  discretization  levels  are  in¬ 
creased  but  instead  let  the  infinite-dimensional  analysis  motivate  potential  sources  of  ill-posedness  in 
the  discrete  least  squares  formulations.  The  effects  of  ill-posedness  are  demonstrated  in  the  context 
of  material  characterization  using  experimental  data. 


•  Nj.  Additionally,  the  iW  x  1 

(22) 
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5  Material  Characterization 


To  illustrate  attributes  of  the  least  squares  parameter  estimation  formulations  (23)  and  (24)  for 
estimating  the  TV  =  TVj  •  Nj  parameters  {v(ECi,  Ejj )},  we  consider  the  characterization  of  the  ferro¬ 
electric  compound  PZT5H  using  experimental  data  which  includes  both  an  unbiased  major  loop  and 
multiple  biased  minor  loops.  Because  data  was  collected  at  0.2  Hz,  relaxation  effects  were  negligible 
so  we  employed  the  piecewise  linear  kernel  (3)  or  (4)  in  the  discretized  model  (10). 

The  unregularized  model  fits  obtained  with  the  discretization  limits  Nt  =  Nj  =  24  (TV  =  576) 
and  TVj  =  Nj  =  48  (TV  =  2304)  using  data  from  all  seven  hysteresis  loops  are  plotted  in  Figure  2 
whereas  those  obtained  using  the  same  quadrature  limits  in  the  regularized  functional  are  given  in 
Figure  3.  Without  regularization,  the  ill-posedness  associated  with  inversion  of  the  discretized  com¬ 
pact  operator  yields  increasingly  inaccurate  model  predictions  as  quadrature  limits  are  increased. 
Regularization  through  the  inclusion  of  the  penalty  term  ([  ]](/][ 2  stabilizes  the  pseudoinverse  by  shift¬ 
ing  singular  values  away  from  the  origin  thus  yielding  the  highly  accurate  fit  observed  in  Figure  3. 
The  regularization  parameter  a  =  5  x  1020  used  to  obtain  the  modeled  behavior  in  Figure  3  was 
computed  using  a  variation  of  the  unbiased  predictive  risk  estimator  (UPRE)  method  discussed  in 
Vogel  [27] .  Smaller  and  larger  values  of  a  yielded  larger  residuals  and  increasingly  inaccurate  modeled 
behavior  —  e.g.,  the  fit  in  Figure  2  corresponds  to  a  =  0. 

Additional  examples  illustrating  the  performance  of  the  model  with  a  priori  functional  choices 
for  and  can  be  found  in  [20]  where  it  is  demonstrated  that  general  densities  are  required  to 
obtain  accurate  characterization  through  the  full  hierarchy  of  drive  conditions. 

6  Concluding  Remarks 

This  paper  addresses  theoretical  and  implementation  issues  associated  with  the  construction  of  mod¬ 
els  used  to  characterize  hysteresis  and  constitutive  nonlinearities  inherent  to  ferroelectric,  ferromag¬ 
netic  and  ferroelastic  materials.  This  unified  characterization  framework  employs  energy  analysis 
to  construct  polarization,  magnetization  and  strain  kernels  at  the  lattice  level.  To  incorporate  the 
effects  of  material  nonhomogeneities,  polycrystallinity,  and  variable  effective  fields,  certain  parame¬ 
ters  are  assumed  to  be  manifestations  of  underlying  distributions  rather  than  constants.  Stochastic 
homogenization  in  this  manner  yields  integral  formulations  having  known  kernels  and  unknown  den- 


Electric  Field  (MV/m)  Electric  Field  (MV/m) 


(a) 


(b) 


Figure  2:  PZT5H  data  and  model  fit  with  general  product  density  v  estimated  using  the  unregularized 
functional  (23)  with  data  from  all  7  loops,  (a)  TV*  =  Nj  =  24  (TV  =  576),  and  (b)  TVj  =  Nj  =  48 
(TV  =  2304). 


11 


Electric  Field  (MV/m)  Electric  Field  (MV/m) 


(a) 


(b) 


Figure  3:  PZT5H  data  and  model  fit  with  general  product  density  v  estimated  using  the  regularized 
Tikhonov  functional  (24)  with  data  from  all  7  loops,  (a)  TVj  =  Nj  =  24  (N  =  576),  and  (b) 
Ni  =  Nj  =  48  (N  =  2304). 


sities  which  must  be  estimated  for  a  given  material  of  transducer  configuration.  Estimation  of  either 
these  general  densities,  or  parameters  in  assumed  functional  forms  for  the  densities,  comprises  the 
parameter  estimation  problem  under  consideration  in  this  paper. 

It  is  noted  that  in  this  homogenized  energy  framework,  one  has  the  option  of  estimating  general 
constituent  densities  u\  and  the  joint  density  u  =  v\  ■  z^2,  or  functional  forms  of  and  V2 
having  a  small  number  of  parameters.  It  is  illustrated  in  [20]  that  the  third  option  leads  to  a  small 
number  of  parameter  to  be  identified  (e.g.,  5  to  6)  but  yields  models  having  limited  accuracy  for 
general  operating  regimes.  Consideration  of  the  joint  density  v  increases  the  dimensionality  of  the 
parameter  estimation  problem  but  yields  a  linear  parameterization  as  required  for  linear  adaptive 
identification  and  control  techniques.  For  this  reason,  we  focused  on  this  formulation  throughout  the 
development. 

It  is  demonstrated  that  the  integral  operators  associated  with  the  model  are  compact  and  hence 
the  inverse  problem  is  ill-posed.  This  motivates  consideration  of  Tikhonov  regularization  to  stabilize 
the  pseudoinverse  by  shifting  singular  values  away  from  the  origin.  It  is  shown  in  the  context  of 
characterizing  PZT5H  material  behavior  that  increasingly  inaccurate  model  fits  are  obtained  as 
discretization  levels  are  increased  in  the  absence  of  regularization  whereas  the  deleterious  effects  of 
ill-posedness  are  avoided  when  suitable  regularization  is  incorporated. 

From  a  practical  perspective,  the  decrease  in  regularity  and  accuracy  due  to  ill-posedness  can  be 
mitigated  by  employing  small  to  moderate  discretization  limits  —  this  forms  the  basis  of  regulariza¬ 
tion  by  coarsening.  Whereas  convergence  issues  are  not  addressed  by  this  latter  approach,  it  yields 
models  which  can  facilitate  control  design. 
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